### r code for the "Birds of a feather" project

library(stm)
library(ggplot2)
library(geomtextpath)
library(ggpubr)
library(ggcorrplot)
library(gridExtra)
library(grid)
library(wordcloud)
library(ggwordcloud)
library(RColorBrewer)
library(tidyverse)
library(tidytext)
library(reshape2)
library(kableExtra)


#### figure 1 #####
## load "hk.xlsx"(sheet name: barometer1)
a<- ggplot(barometer1, aes(x=year, y=percent, group = type)) + geom_line(aes(color = type, linetype = type),linewidth=0.8) + geom_point(aes(shape=type),size=2) + 
  xlab("") + ylab("Percentage") + scale_x_discrete(guide= guide_axis(n.dodge = 2)) +
  theme(text = element_text(size = 12),
        legend.position = c(.8,.9), legend.title = element_blank(),
        strip.text.x = element_text(size = 12),
        axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12),
        plot.title = element_text(hjust = 0.5, size = 12, face = "plain"),
        panel.grid = element_blank(), 
        panel.background = element_rect(color = 'black', fill = 'transparent')) +
  ggtitle("Economy vs Democracy")

## load "hk.xlsx"(sheet name: barometer2)
barometer2 <- transform(barometer2, age_group = reorder(age_group, number))
b<- ggplot(barometer2, aes(x=year, y=percent, group = type)) + geom_line(aes(color = type, linetype = type),linewidth=0.8) + geom_point(aes(shape=type),size=2) + 
  xlab("") + ylab("Percentage") + scale_x_discrete(guide= guide_axis(n.dodge = 2)) + facet_grid(age_group ~ .) +
  theme(text = element_text(size = 12),
        legend.position = "none", legend.title = element_blank(),
        strip.text.y = element_text(size = 15),
        axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12),
        plot.title = element_text(hjust = 0.5, size = 12, face = "plain"),
        panel.grid = element_blank(), 
        panel.background = element_rect(color = 'black', fill = 'transparent')) +
  ggtitle("Economy vs Democracy (across Age Groups)")

grid.arrange(a,b,ncol = 2)



#### figure 2 #####
## load "words.xlsx" (sheet name: tf)
a<-ggplot(words, aes(label = visi, size = tf_idf_v)) + geom_text_wordcloud_area() + scale_size_area(max_size = 15) + theme_minimal()

words <- transform(words, words_v2 = reorder(words_v2, number))
b<-ggplot(words, aes(x=words_v2, y=tf_idf_v)) + geom_bar(stat='identity', size=0.5, fill = "blue") + xlab("") + ylab("TF-IDF Score") + 
  geom_text(aes(label=words_v2), hjust=-.1, size = 2.5) + 
  theme(axis.text.y = element_blank(),
        axis.ticks.y=element_blank(),
        plot.caption = element_text(hjust = 0),
        panel.grid = element_blank(), 
        panel.background = element_rect(color = 'black', fill = 'transparent')) + coord_flip() +ylim(c(0,0.1))

grid.arrange(a,b, ncol = 2)

#### figure 3 #####
## load "hk_migrant_fenci.csv"
processed <- textProcessor(hk_migrant$content, metadata = hk_migrant, lowercase = FALSE, removestopwords = FALSE, 
                           removenumbers = FALSE, removepunctuation = TRUE, ucp = FALSE, stem = FALSE, sparselevel = 1, language = "na", 
                           wordLengths = c(2, Inf))

plotRemoved(processed$documents, lower.thresh=seq(1,200,by=10))

out <- prepDocuments(processed$documents, processed$vocab, 
                     processed$meta, lower.thresh = 15)

docs <- out$documents
vocab <- out$vocab
meta <-out$meta

select <- searchK(documents = out$documents, vocab = out$vocab, K = c(3, 4, 5, 6, 7, 8, 9, 10), 
                  max.em.its = 25, prevalence = ~ newspaper + s(date), data = out$meta)

set.seed(12345)
hk_migrant_stm <- stm(documents = out$documents, vocab = out$vocab, K = 5, prevalence = ~ newspaper + s(date), max.em.its = 25, data = out$meta,
                      init.type = "Spectral", verbose = TRUE)

labelTopics(hk_migrant_stm, c(1:5), n = 10, frexweight = 0.5)

plot(hk_migrant_stm, type="summary", xlim=c(0, 0.5), main = "", 
     custom.labels = c("Social Welfare", 
                       "Hong Kong-Mainland Economic & Political Tensions", 
                       "Education", 
                       "Housing Problems", 
                       "Hong Kong Elections"), cex.lab = 1.0, text.cex = 1.0)

thoughts1 <- findThoughts(hk_migrant_stm, texts = out$meta$title, n = 30, topics = 1)$docs[[1]]


#### figure 4 #####
prep <- estimateEffect(~ newspaper + s(date), hk_migrant_stm, meta = out$meta, uncertainty = "Global")
par(mfrow=c(3,2))

plot(prep, "date", method = "continuous", topics = c(1), printlegend = F, xaxt = "n", ci.level = T, linecol = "black", 
     ylab = "", xlab = "", ylim = c(0, 0.5), axes=TRUE)
title("Topic 1: Social Welfare", line = -1, cex.main = 1, font.main= 1)
axis(1, at = c(20030201, 20070101, 20110101, 20150101, 20190101), 
     labels=c("2003-02", "2007-01", "2011-01", "2015-01","2019-01"), font=0.1, cex= 0.3)

plot(prep, "date", method = "continuous", topics = c(2), printlegend = F, xaxt = "n", ci.level = T, linecol = "black", 
     ylab = "", xlab = "", ylim = c(0, 0.5), axes=TRUE)
title("Topic 2: HK-Mainland Tensions", line = -1, cex.main = 1, font.main= 1)
axis(1, at = c(20030201, 20070101, 20110101, 20150101, 20190101), 
     labels=c("2003-02", "2007-01", "2011-01", "2015-01","2019-01"), font=0.1, cex= 0.3)

plot(prep, "date", method = "continuous", topics = c(3), printlegend = F, xaxt = "n", ci.level = T, linecol = "black", 
     ylab = "", xlab = "", ylim = c(0, 0.5), axes=TRUE)
title("Topic 3: Education", line = -1, cex.main = 1, font.main= 1)
axis(1, at = c(20030201, 20070101, 20110101, 20150101, 20190101), 
     labels=c("2003-02", "2007-01", "2011-01", "2015-01","2019-01"), font=0.1, cex= 0.3)

plot(prep, "date", method = "continuous", topics = c(4), printlegend = F, xaxt = "n", ci.level = T, linecol = "black", 
     ylab = "", xlab = "", ylim = c(0, 0.5), axes=TRUE)
title("Topic 4: Housing Problems", line = -1, cex.main = 1, font.main= 1)
axis(1, at = c(20030201, 20070101, 20110101, 20150101, 20190101), 
     labels=c("2003-02", "2007-01", "2011-01", "2015-01","2019-01"), font=0.1, cex= 0.3)

plot(prep, "date", method = "continuous", topics = c(5), printlegend = F, xaxt = "n", ci.level = T, linecol = "black", 
     ylab = "", xlab = "", ylim = c(0, 0.5), axes=TRUE)
title("Topic 5: Hong Kong Elections", line = -1, cex.main = 1, font.main= 1)
axis(1, at = c(20030201, 20070101, 20110101, 20150101, 20190101), 
     labels=c("2003-02", "2007-01", "2011-01", "2015-01","2019-01"), font=0.1, cex= 0.3)


#### figure 6 #####
## t test results from Stata
## load "hk.xlsx"(sheet name: attitude_ttest)
theme_set(theme_classic())
attitude_ttest <- transform(attitude_ttest, group = reorder(group, number))
p<- ggplot(attitude_ttest, aes(x=group, y=means, fill = group)) + scale_fill_brewer(palette="Reds") +
  geom_bar(aes(x=group, y=means), stat="identity", color="black" , size = .5) + 
  geom_errorbar(aes(ymin=means-sd, ymax=means+sd), size = .5, width=.2) + xlab("") + ylab("")+ ylim(0,5)+
  annotate(geom = "path", x=c(1,1,2,2),y=c(3.5,3.6,3.6,3.5)) +
  annotate(geom="text", x=1.5, y=3.8, label="p = 0.869", size = 3.5, fontface="plain") + 
  annotate(geom = "path", x=c(1,1,3,3),y=c(3.9,4,4,3.7)) +
  annotate(geom="text", x=2, y=4.2, label="p = 0.266", size = 3.5, fontface="plain") + 
  annotate(geom = "path", x=c(1,1,4,4),y=c(4.3,4.4,4.4,3.6)) +
  annotate(geom="text", x=2.5, y=4.6, label="p = 0.541", size = 3.5, fontface="plain") +
  annotate(geom = "path", x=c(1,1,5,5),y=c(4.7,4.8,4.8,3.7)) +
  annotate(geom="text", x=3, y=5, label="p = 0.453", size = 3.5, fontface="plain") +
  theme(legend.position = "none",
        axis.text = element_text(size=12, color="black"),
        panel.grid = element_blank(), 
        panel.background = element_rect(color = 'black', fill = 'transparent'),
        plot.title = element_text(size = 12, face = "plain")) + ggtitle("Panel A: Attitude toward Immigrants")

## load "hk.xlsx"(sheet name: policy_ttest)
policy_ttest <- transform(policy_ttest, group = reorder(group, number))
q<- ggplot(policy_ttest, aes(x=group, y=means, fill = group)) + scale_fill_brewer(palette="Blues") +
  geom_bar(stat="identity", color="black" , size = .5) + 
  geom_errorbar(aes(ymin=means-sd, ymax=means+sd), size = .5, width=.2) + xlab("") + ylab("")+ ylim(0,6.5)+
  annotate(geom = "path", x=c(1,1,2,2),y=c(4.8,4.9,4.9,4.8)) +
  annotate(geom="text", x=1.5, y=5.1, label="p = 0.698", size = 3.5, fontface="plain") + 
  annotate(geom = "path", x=c(1,1,3,3),y=c(5.2,5.3,5.3,4.8)) +
  annotate(geom="text", x=2, y=5.5, label="p = 0.965", size = 3.5, fontface="plain") +
  annotate(geom = "path", x=c(1,1,4,4),y=c(5.6,5.7,5.7,4.8)) +
  annotate(geom="text", x=2.5, y=5.9, label="p = 0.736", size = 3.5, fontface="plain") +
  annotate(geom = "path", x=c(1,1,5,5),y=c(6,6.1,6.1,4.8), colour = 'red') +
  annotate(geom="text", x=3, y=6.3, label="p = 0.058", size = 3.5, fontface="plain", colour = 'red') +
  theme(legend.position = "none",
        axis.text = element_text(size=12, color="black"),
        plot.caption = element_text(hjust = 0),
        panel.grid = element_blank(), 
        panel.background = element_rect(color = 'black', fill = 'transparent'),
        plot.title = element_text(size = 12, face = "plain")) + ggtitle("Panel B: Support Restriction") 

grid.arrange(p,q)


#### figure 7 ####
## load "hk.xlsx"(sheet name: main_reg)
main_reg <- transform(main_reg, type=reorder(type, number))
main_reg <- transform(main_reg, group=reorder(group, number2))
ggplot(main_reg, aes(x=coef, y=type)) + geom_point(size=4)+ geom_errorbar(aes(xmin=low90, xmax=high90), width=0, size=2, color="blue4") + 
  geom_errorbar(aes(xmin=low95, xmax=high95), size =1, width=0, color="blue") +
  ylab("") + xlab("Regression Result") + geom_vline(xintercept = 0, linetype="dotted", color = "red") +
  facet_wrap(~group, scales = "free")+
  theme(panel.grid = element_blank(), 
        panel.background = element_rect(color = 'black', fill = 'transparent'),
        strip.background = element_rect(colour="black",fill="white"),
        axis.text = element_text(size = 12),
        strip.text.x = element_text(size = 12))

#### figure 8 ####
## subgroup t test result from Stata
## load "hk.xlsx"(sheet name: group_politics)
group_politics <- transform(group_politics, group = reorder(group, number))

text5_income_1 <- data.frame(label = c("p = 0.844", "p = 0.045", "p = 0.506"),
                             income = c("1 low", "2 middle", "3 high"),
                             x     = c(1.5, 1.5, 1.5),
                             y     = c(5.1, 5.1, 5.1))
group_politics_ <- group_politics[rowSums(is.na(group_politics)) == 0,] 
a<- ggbarplot(group_politics_, x = "group", y = "policy", fill = "group", palette="Blues",
              add = "mean_se", facet.by = "income", short.panel.labs = FALSE)+ 
  geom_text(data = text5_income_1, mapping = aes(x = x, y = y, label = label))+
  theme(legend.position = "none",
        strip.text.x = element_text(size = 12),
        axis.text = element_text(size=12, color="black"),
        plot.title = element_text(size = 12, face = "plain"),
        axis.title.y = element_text(size = 12)) + xlab("") + ylab("Support Restriction") + ggtitle("Panel A: By Income")



text5_born_1 <- data.frame(label = c("p = 0.036", "p = 0.995", "p = 0.840"),
                           born = c("1 Hong Kong", "2 Mainland China", "3 Overseas"),
                           x     = c(1.5, 1.5, 1.5),
                           y     = c(5.5, 5.5, 6.2))

b<- ggbarplot(group_politics, x = "group", y = "policy", fill = "group", palette="Blues",
              add = "mean_se", facet.by = "born", short.panel.labs = FALSE)+ 
  geom_text(data = text5_born_1, mapping = aes(x = x, y = y, label = label))+
  theme(legend.position = "none",
        strip.text.x = element_text(size = 12),
        axis.text = element_text(size=12, color="black"),
        plot.title = element_text(size = 12, face = "plain"),
        axis.title.y = element_text(size = 12)) + xlab("") + ylab("Support Restriction") + ggtitle("Panel B: By Birthplace")


text5_edu_1 <- data.frame(label = c("p = 0.335", "p = 0.091"),
                          edu = c("1 non-college", "2 college or above"),
                          x = c(1.5, 1.5),
                          y = c(5.2, 5.2))

c <- ggbarplot(group_politics, x = "group", y = "policy", fill = "group", palette="Blues",
               add = "mean_se", facet.by = "edu", short.panel.labs = FALSE)+ 
  geom_text(data = text5_edu_1, mapping = aes(x = x, y = y, label = label))+
  theme(legend.position = "none",
        strip.text.x = element_text(size = 12),
        axis.text = element_text(size=12, color="black"),
        plot.title = element_text(size = 12, face = "plain"),
        axis.title.y = element_text(size = 12)) + xlab("") + ylab("Support Restriction") + ggtitle("Panel C: By Education")


grid.arrange(a,b,c)


################### appendix #########################

#### figure A1 ####
## load "hk.xlsx"(sheet name: mainland_immigrant)
a<- ggplot() +
  geom_line(mapping = aes(x = mainland_immigrant$year, y = mainland_immigrant$number), size = 1.5, color = "steelblue") +
  geom_text(aes(x = mainland_immigrant$year, y = mainland_immigrant$number, label=mainland_immigrant$number), vjust=1.6, color="black", size=4) +
  theme_bw() + xlab("Year") + ylab("Number of Mainland Immigrants") +
  scale_x_continuous(breaks = c(1997,1999,2001,2003,2005,2007,2009,2011,2013,2015,2017,2019))

b<- ggplot() +
  geom_bar(mapping = aes(x = mainland_immigrant$year, y = mainland_immigrant$ratio), stat="identity", fill="steelblue") +
  geom_text(aes(x = mainland_immigrant$year, y = mainland_immigrant$ratio, label=sprintf("%.2f", mainland_immigrant$ratio)), vjust=1.6, color="black", size=4) +
  theme_bw() + xlab("Year") + ylab("Ratio") +
  scale_x_continuous(breaks = c(1997,1999,2001,2003,2005,2007,2009,2011,2013,2015,2017,2019))

grid.arrange(a,b, ncol = 2)


#### figure A2 ####
## load "hk.xlsx"(sheet name: hk_identity)
mytype <- hk_identity$Type

ggplot(data=hk_identity, aes(x=Year, y=Ratio, group=mytype)) + geom_line(aes(linetype=mytype, color=mytype), size=1.0) + theme_bw() + 
  theme(legend.position = c(0.15, 0.85), legend.title = element_blank()) + ylab("Percentage") 


#### figure A3 ####
ggplot(hk_attitude, aes(x=year, y=government, group=type1)) + geom_line(aes(color=type1, linetype=type1), size=0.8) + xlab("Year") + ylab("Percentage") +
  theme(legend.position = "top", legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 10, face = "plain"),
        panel.grid = element_blank(), 
        panel.background = element_rect(color = 'black', fill = 'transparent')) +
  ggtitle("Mainland Government")

ggplot(hk_attitude2, aes(x=year2, y=people, group=type2)) + geom_line(aes(color=type2, linetype=type2), size=0.8) + xlab("Year") + ylab("Percentage") +
  theme(legend.position = "top", legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 10, face = "plain"),
        panel.grid = element_blank(), 
        panel.background = element_rect(color = 'black', fill = 'transparent')) +
  ggtitle("Mainland Chinese")


##### figure B1 #####
## load "hk.xlsx"(sheet name: news_report)
theme_set(theme_bw())
ggplot(news_report, aes(x= year, y=number, color=type, lty =type)) + geom_line(size = 0.8) + xlab("") + ylab("Number of News Articles") + 
  theme(legend.position = c(0.2, 0.9), legend.title=element_blank())



#### figure B2 ######
## load "hk.xlsx"(sheet name: top_model)
ggplot(top_model, aes(x=semantic, exclusivity)) + geom_point() + geom_text(aes(label=topic_number),hjust=0.5, vjust=-0.5) +
  xlab("Semantic Coherence") + ylab("Exclusivity") +
  theme(panel.grid = element_blank(), 
        panel.background = element_rect(color = 'black', fill = 'transparent')) +
  #add an arrow
  geom_segment(aes(x = -53.88869, y = 8.6, xend = -53.88869, yend = 8.8), arrow = arrow(length = unit(0.3, "cm")), colour = "red")


#### figure B3 is plotted using Stata


#### figure B4-B6 ####
## load "hk.xlsx"(sheet name: hk_negative)
Sys.setlocale("LC_TIME", "C")

hk_negative$sdate <- as.Date(hk_negative$sdate, format = "%Y-%m-%d")
ggplot(hk_negative, aes(x = sdate, y = negative_weighted)) + 
  geom_line(size = 0.5) + scale_x_date(date_labels="%Y-%m",date_breaks  ="36 month") +
  ggtitle("Negative Sentiment Dynamics") + xlab("") + ylab("") +
  theme_bw()+
  theme(plot.title = element_text(hjust = 0.5, size = 10, face = "plain"), 
        legend.title = element_text(size=8), 
        legend.text = element_text(size=8),
        legend.position = "none") +
  facet_wrap(~topic, ncol = 3)

# fear distribution
## load "hk.xlsx"(sheet name: hk_fear)
hk_fear$sdate <- as.Date(hk_fear$sdate , format = "%Y-%m-%d")
ggplot(hk_fear, aes(x = sdate, y = fear_weighted)) + 
  geom_line(size = 0.5) + scale_x_date(date_labels="%Y-%m",date_breaks  ="36 month") +
  ggtitle("Fear Dynamics") + xlab("") + ylab("") +
  theme_bw()+
  theme(plot.title = element_text(hjust = 0.5, size = 10, face = "plain"), 
        legend.title = element_text(size=8), 
        legend.text = element_text(size=8),
        legend.position = "none") +
  facet_wrap(~topic, ncol = 3)

# disgust distribution
## load "hk.xlsx"(sheet name: hk_disgust)
hk_disgust$sdate <- as.Date(hk_disgust$sdate, format = "%Y-%m-%d")
ggplot(hk_disgust, aes(x = sdate, y = disgust_weighted)) + 
  geom_line(size = 0.5) + scale_x_date(date_labels="%Y-%m",date_breaks  ="36 month") +
  ggtitle("Disgust Dynamics") + xlab("") + ylab("") +
  theme_bw()+
  theme(plot.title = element_text(hjust = 0.5, size = 10, face = "plain"), 
        legend.title = element_text(size=8), 
        legend.text = element_text(size=8),
        legend.position = "none") +
  facet_wrap(~topic, ncol = 3)


#### figure B7 ####
## load "words.xlsx" (sheet name: tf)
c<- ggplot(words, aes(x=frequency_v, y=words_v1)) + geom_bar(stat='identity', size=0.5, fill = "blue") + xlab("words frequency") + ylab("") +
  ggtitle("news articles on Mainland visitors") 

d<- ggplot(words, aes(x=tf_idf_v, y=words_v2)) + geom_bar(stat='identity', size=0.5, fill = "blue") + xlab("TF-IDF Score") + ylab("") +
  ggtitle("news articles on Mainland visitors")

grid.arrange(c,d, ncol = 2)


#### figure B8 ####
## load "hk_visitor_fenci.csv"
processed_v <- textProcessor(hk_visitor$content, metadata = hk_visitor, 
                             lowercase = FALSE, removestopwords = FALSE,
                             removenumbers = FALSE, removepunctuation = TRUE,
                             ucp = FALSE, stem = FALSE, 
                             sparselevel = 1, language = "na", 
                             wordLengths = c(2, Inf))

plotRemoved(processed_v$documents, lower.thresh=seq(1,200,by=10))

out_v <- prepDocuments(processed_v$documents, processed_v$vocab, 
                       processed_v$meta, lower.thresh = 15)

docs_v <- out_v$documents
vocab_v <- out_v$vocab
meta_v <-out_v$meta

hk_visitor_stm <- stm(documents = out_v$documents, vocab = out_v$vocab,
                      K = 4, 
                      max.em.its = 25, data = out_v$meta,
                      init.type = "Spectral", verbose = TRUE)

gamma <- tidy(hk_visitor_stm, matrix = "gamma", document_names = rownames(hk_visitor))

gamma <- gamma %>%
  group_by(topic) %>%
  summarise(gamma = mean(gamma)) %>%
  arrange(desc(gamma)) %>%
  
  mutate(topic = paste0("Topic ", topic),
         topic = reorder(topic, gamma))

gamma %>%
  select(topic, gamma) %>%
  kable(digits = 3, 
        col.names = c("Topic", "Expected topic proportion"))
rm(keywards)
proportion <- c(0.322, 0.280, 0.216, 0.182)
name <- c("Topic 1: \n Mainlanders \n Traveling to HK", "Topic 3: \n Official Statement", 
          "Topic 2: \n Anti-ELAB Movement", "Topic 4: \n Panic Buying \n Milk Powder")
keywords <- c("Entry   Shopping city   Port", "OCTS   Autonomy   Constitution", "Rioter   Demonstrators   Government House",
              "Milk Powder   Milk Powder Restriction   Supplying")
number <- c(4,3,2,1)
lvyou <- data.frame(proportion, name, keywords, number)

lvyou <- transform(lvyou, name = reorder(name, number))
ggplot(lvyou, aes(x = proportion, y = name)) + geom_segment(aes(x=0, xend=proportion, y=name, yend=name)) +
  geom_text(aes(label = keywords), size = 8, hjust = -0.05, vjust = 0.5)+ xlim(0,0.5) + ylab("")+ 
  theme(panel.background = element_rect(fill = "white", colour = "black"), axis.text = element_text(size = 15))



##### figure C1 ####
## load "hk.xlsx"(sheet name: interactive)
ggplot(interactive, aes(x= factor(experience, level = c('0 (very bad)', '1', '2', '3', '4', '5', '6', '7', '8', '9', '10 (very good)')), 
               y=coef, group=1)) + geom_point(size = 2) + geom_line() + 
  geom_ribbon(aes(ymin=low.ci, ymax=high.ci), linetype=2, alpha=0.1) +
  geom_hline(yintercept=0, linetype="dashed", color = "red") +
  labs(x="Interactive Experience with Mainlanders", y="Support Restriction (Linear Prediction)") +
  theme(plot.caption = element_text(hjust = 0),
        axis.text = element_text(size=10, color="black"),
        panel.grid = element_blank(), 
        panel.background = element_rect(color = 'black', fill = 'transparent'))


#### figure C2-C3 #####
## t test results from Stata

### education
## load "hk.xlsx"(sheet name: group_econ)
text2_edu2 <- data.frame(label = c("p = 0.521", "p = 0.829"),
                         education = c("1 below bachelor", "2 bachelor or above"),
                         x     = c(1.5, 1.5),
                         y     = c(4.1, 4.1))

a1<- ggbarplot(group_econ, x = "group", y = "attitude", fill = "group", palette="Blues", 
               add = "mean_se", facet.by = "education", short.panel.labs = FALSE)+ 
  geom_text(data = text2_edu2, mapping = aes(x = x, y = y, label = label))+
  theme(legend.position = "none",
        axis.text = element_text(size=10, color="black"),
        plot.title = element_text(size = 10, face = "plain"),
        axis.title.y = element_text(size = 10)) + ggtitle("Panel A: Economic Contribution") +
  xlab("") + ylab("Attitude")


text2_edu <- data.frame(label = c("p = 0.054", "p = 0.129"),
                        education = c("1 below bachelor", "2 bachelor or above"),
                        x     = c(1.5, 1.5),
                        y     = c(5.5, 5.5))

a2<- ggbarplot(group_econ, x = "group", y = "policy", fill = "group", palette="Blues", 
               add = "mean_se", facet.by = "education", short.panel.labs = FALSE)+ 
  geom_text(data = text2_edu, mapping = aes(x = x, y = y, label = label))+
  theme(legend.position = "none",
        axis.text = element_text(size=10, color="black"),
        plot.title = element_text(size = 10, face = "plain"),
        axis.title.y = element_text(size = 10)) + ggtitle("Panel A: Economic Contribution") +
  xlab("") + ylab("Support Restriction")


### education
## load "hk.xlsx"(sheet name: group_social)
text3_edu <- data.frame(label = c("p = 0.166", "p = 0.098"),
                        education = c("1 below bachelor", "2 bachelor or above"),
                        x     = c(1.5, 1.5),
                        y     = c(4.2, 4.2))

group_social$group <- factor(group_social$group, levels = c("control group", "social norms"))

b1<- ggbarplot(group_social, x = "group", y = "attitude", fill = "group", palette="Blues",
               add = "mean_se", facet.by = "education", short.panel.labs = FALSE)+ 
  geom_text(data = text3_edu, mapping = aes(x = x, y = y, label = label))+
  theme(legend.position = "none",
        axis.text = element_text(size=10, color="black"),
        plot.title = element_text(size = 10, face = "plain"),
        axis.title.y = element_text(size = 10)) + ggtitle("Panel B: Social Norms") +
  xlab("") + ylab("Attitude")


text3_edu2 <- data.frame(label = c("p = 0.171", "p = 0.553"),
                         education = c("1 below bachelor", "2 bachelor or above"),
                         x     = c(1.5, 1.5),
                         y     = c(5.4, 5.2))

b2<- ggbarplot(group_social, x = "group", y = "policy", fill = "group", palette="Blues",
               add = "mean_se", facet.by = "education", short.panel.labs = FALSE)+ 
  geom_text(data = text3_edu2, mapping = aes(x = x, y = y, label = label))+
  theme(legend.position = "none",
        axis.text = element_text(size=10, color="black"),
        plot.title = element_text(size = 10, face = "plain"),
        axis.title.y = element_text(size = 10)) + ggtitle("Panel B: Social Norms") +
  xlab("") + ylab("Support Restriction")


### income
## load "hk.xlsx"(sheet name: group_family)
text4_income <- data.frame(label = c("p = 0.313", "p = 0.879", "p = 0.065"),
                           income = c("1 low", "2 middle", "3 high"),
                           x     = c(1.5, 1.5, 1.5),
                           y     = c(4.2, 4.2, 4.2))

group_family$group <- factor(group_family$group, levels = c("control group", "family history"))

c1<- ggbarplot(group_family, x = 'group', y = "attitude", fill = "group", palette="Blues", 
               add = "mean_se", facet.by = "income", short.panel.labs = FALSE)+ 
  geom_text(data = text4_income, mapping = aes(x = x, y = y, label = label))+
  theme(legend.position = "none",
        axis.text = element_text(size=10, color="black"),
        plot.title = element_text(size = 10, face = "plain"),
        axis.title.y = element_text(size = 10)) + ggtitle("Panel C: Family History") +
  xlab("") + ylab("Attitude")


text4_income2 <- data.frame(label = c("p = 0.737", "p = 0.937", "p = 0.516"),
                            income = c("1 low", "2 middle", "3 high"),
                            x     = c(1.5, 1.5, 1.5),
                            y     = c(5.2, 5.2, 5.4))

c2<- ggbarplot(group_family, x = 'group', y = "policy", fill = "group", palette="Blues",
               add = "mean_se", facet.by = "income", short.panel.labs = FALSE)+ 
  geom_text(data = text4_income2, mapping = aes(x = x, y = y, label = label))+
  theme(legend.position = "none",
        axis.text = element_text(size=10, color="black"),
        plot.title = element_text(size = 10, face = "plain"),
        axis.title.y = element_text(size = 10)) + ggtitle("Panel C: Family History") +
  xlab("") + ylab("Support Restriction") 


grid.arrange(a1,b1,c1)

grid.arrange(a2,b2,c2)




